!===============================================================================
! Bulk Aerosol Model
!===============================================================================
module aero_model
  use shr_kind_mod,      only: r8 => shr_kind_r8
  use constituents,      only: pcnst, cnst_name, cnst_get_ind
  use ppgrid,            only: pcols, pver, pverp
  use cam_abortutils,    only: endrun
  use cam_logfile,       only: iulog
  use perf_mod,          only: t_startf, t_stopf
  use camsrfexch,        only: cam_in_t, cam_out_t
  use aerodep_flx,       only: aerodep_flx_prescribed
  use physics_types,     only: physics_state, physics_ptend, physics_ptend_init
  use physics_buffer,    only: physics_buffer_desc
  use physconst,         only: gravit, rair
  use dust_model,        only: dust_active, dust_names, dust_nbin
  use seasalt_model,     only: sslt_active=>seasalt_active, seasalt_names, seasalt_nbin
  use spmd_utils,        only: masterproc
  use physics_buffer,    only: pbuf_get_field, pbuf_get_index
  use cam_history,       only: outfld
  use infnan,            only: nan, assignment(=)

  implicit none
  private

  public :: aero_model_readnl
  public :: aero_model_register
  public :: aero_model_init
  public :: aero_model_gasaerexch ! create, grow, change, and shrink aerosols.
  public :: aero_model_drydep     ! aerosol dry deposition and sediment
  public :: aero_model_wetdep     ! aerosol wet removal
  public :: aero_model_emissions  ! aerosol emissions
  public :: aero_model_surfarea    ! tropospheric aerosol wet surface area for chemistry
  public :: aero_model_strat_surfarea   ! stub

 ! Misc private data

  integer :: so4_ndx, cb2_ndx, oc2_ndx, nit_ndx
  integer :: soa_ndx, soai_ndx, soam_ndx, soab_ndx, soat_ndx, soax_ndx

  ! Namelist variables
  character(len=16) :: wetdep_list(PCNST) = ' '
  character(len=16) :: drydep_list(PCNST) = ' '

  integer :: ndrydep = 0
  integer, allocatable :: drydep_indices(:)
  integer :: nwetdep = 0
  integer, allocatable :: wetdep_indices(:)
  logical :: drydep_lq(PCNST)
  logical :: wetdep_lq(PCNST)

  integer :: fracis_idx = 0

  real(r8) aer_sol_facti(PCNST) ! in-cloud solubility factor
  real(r8) aer_sol_factb(PCNST) ! below-cloud solubility factor
  real(r8) aer_scav_coef(PCNST)

contains

  !=============================================================================
  ! reads aerosol namelist options
  !=============================================================================
  subroutine aero_model_readnl(nlfile)

    use namelist_utils,  only: find_group_name
    use units,           only: getunit, freeunit
    use mpishorthand

    character(len=*), intent(in) :: nlfile  ! filepath for file containing namelist input

    ! Local variables
    integer :: unitn, ierr
    character(len=*), parameter :: subname = 'aero_model_readnl'

    ! Namelist variables
    character(len=16) :: aer_wetdep_list(PCNST) = ' '
    character(len=16) :: aer_drydep_list(PCNST) = ' '

    namelist /aerosol_nl/ aer_wetdep_list, aer_drydep_list
    namelist /aerosol_nl/ aer_sol_facti, aer_sol_factb, aer_scav_coef

    aer_sol_facti = nan
    aer_sol_factb = nan
    aer_scav_coef = nan

    ! Read namelist
    if (masterproc) then
       unitn = getunit()
       open( unitn, file=trim(nlfile), status='old' )
       call find_group_name(unitn, 'aerosol_nl', status=ierr)
       if (ierr == 0) then
          read(unitn, aerosol_nl, iostat=ierr)
          if (ierr /= 0) then
             call endrun(subname // ':: ERROR reading namelist')
          end if
       end if
       close(unitn)
       call freeunit(unitn)
    end if

#ifdef SPMD
    ! Broadcast namelist variables
    call mpibcast(aer_wetdep_list, len(aer_wetdep_list(1))*pcnst, mpichar, 0, mpicom)
    call mpibcast(aer_drydep_list, len(aer_drydep_list(1))*pcnst, mpichar, 0, mpicom)
    call mpibcast(aer_sol_facti, pcnst, mpir8, 0, mpicom)
    call mpibcast(aer_sol_factb, pcnst, mpir8, 0, mpicom)
    call mpibcast(aer_scav_coef, pcnst, mpir8, 0, mpicom)
#endif

    wetdep_list = aer_wetdep_list
    drydep_list = aer_drydep_list

  end subroutine aero_model_readnl

  !=============================================================================
  !=============================================================================
  subroutine aero_model_register()
    use mo_setsoa, only : soa_register

    call soa_register()
  end subroutine aero_model_register

  !=============================================================================
  !=============================================================================
  subroutine aero_model_init( pbuf2d )

    use mo_chem_utls,  only: get_inv_ndx, get_spc_ndx
    use cam_history,   only: addfld, add_default, horiz_only
    use phys_control,  only: phys_getopts
    use mo_aerosols,   only: aerosols_inti
    use mo_setsoa,     only: soa_inti
    use dust_model,    only: dust_init
    use seasalt_model, only: seasalt_init
    use drydep_mod,    only: inidrydep
    use wetdep,        only: wetdep_init
    use mo_setsox,     only: has_sox

    ! args
    type(physics_buffer_desc), pointer :: pbuf2d(:,:)

    ! local vars
    character(len=12), parameter :: subrname = 'aero_model_init'
    integer :: m, id
    character(len=20) :: dummy
    logical  :: history_aerosol ! Output MAM or SECT aerosol tendencies
    logical  :: history_dust    ! Output dust

    call phys_getopts( history_aerosol_out = history_aerosol,&
                       history_dust_out    = history_dust   )

    call aerosols_inti()
    call soa_inti(pbuf2d)
    call dust_init()
    call seasalt_init()
    call wetdep_init()

    fracis_idx = pbuf_get_index('FRACIS')

    nwetdep = 0
    ndrydep = 0

    count_species: do m = 1,pcnst
       if ( len_trim(wetdep_list(m)) /= 0 ) then
          nwetdep = nwetdep+1
       endif
       if ( len_trim(drydep_list(m)) /= 0 ) then
          ndrydep = ndrydep+1
       endif
    enddo count_species

    if (nwetdep>0) &
         allocate(wetdep_indices(nwetdep))
    if (ndrydep>0) &
         allocate(drydep_indices(ndrydep))

    do m = 1,ndrydep
       call cnst_get_ind ( drydep_list(m), id, abort=.false. )
       if (id>0) then
          drydep_indices(m) = id
       else
          call endrun(subrname//': invalid drydep species: '//trim(drydep_list(m)) )
       endif

       if (masterproc) then
          write(iulog,*) subrname//': '//drydep_list(m)//' will have drydep applied'
       endif
    enddo
    do m = 1,nwetdep
       call cnst_get_ind ( wetdep_list(m), id, abort=.false. )
       if (id>0) then
          wetdep_indices(m) = id
       else
          call endrun(subrname//': invalid wetdep species: '//trim(wetdep_list(m)) )
       endif

       if (masterproc) then
          write(iulog,*) subrname//': '//wetdep_list(m)//' will have wet removal'
       endif
    enddo

    ! set flags for drydep tendencies
    drydep_lq(:) = .false.
    do m=1,ndrydep
       id = drydep_indices(m)
       drydep_lq(id) =  .true.
    enddo

    ! set flags for wetdep tendencies
    wetdep_lq(:) = .false.
    do m=1,nwetdep
       id = wetdep_indices(m)
       wetdep_lq(id) = .true.
    enddo

    do m = 1,ndrydep

       dummy = trim(drydep_list(m)) // 'TB'
       call addfld (dummy,horiz_only, 'A','kg/m2/s',trim(drydep_list(m))//' turbulent dry deposition flux')
       if ( history_aerosol ) then
          call add_default (dummy, 1, ' ')
       endif
       dummy = trim(drydep_list(m))  // 'GV'
       call addfld (dummy,horiz_only, 'A','kg/m2/s',trim(drydep_list(m)) //' gravitational dry deposition flux')
       if ( history_aerosol ) then
          call add_default (dummy, 1, ' ')
       endif
       dummy = trim(drydep_list(m))  // 'DD'
       call addfld (dummy,horiz_only, 'A','kg/m2/s',trim(drydep_list(m)) //' dry deposition flux at bottom (grav + turb)')
       if ( history_aerosol ) then
          call add_default (dummy, 1, ' ')
       endif
       dummy = trim(drydep_list(m)) // 'DT'
       call addfld (dummy,(/ 'lev' /), 'A','kg/kg/s',trim(drydep_list(m))//' dry deposition')
       if ( history_aerosol ) then
          call add_default (dummy, 1, ' ')
       endif
       dummy = trim(drydep_list(m)) // 'DV'
       call addfld (dummy,(/ 'lev' /), 'A','m/s',trim(drydep_list(m))//' deposition velocity')
       if ( history_aerosol ) then
          call add_default (dummy, 1, ' ')
       endif

    enddo

    if (ndrydep>0) then

       call inidrydep(rair, gravit)

       dummy = 'RAM1'
       call addfld (dummy,horiz_only, 'A','frac','RAM1')
       if ( history_aerosol ) then
          call add_default (dummy, 1, ' ')
       endif
       dummy = 'airFV'
       call addfld (dummy,horiz_only, 'A','frac','FV')
       if ( history_aerosol ) then
          call add_default (dummy, 1, ' ')
       endif

       if (sslt_active) then
          dummy = 'SSTSFDRY'
          call addfld (dummy,horiz_only, 'A','kg/m2/s','Sea salt deposition flux at surface')
          if ( history_aerosol ) then
             call add_default (dummy, 1, ' ')
          endif
       endif
       if (dust_active) then
          dummy = 'DSTSFDRY'
          call addfld (dummy,horiz_only, 'A','kg/m2/s','Dust deposition flux at surface')
          if ( history_aerosol ) then
             call add_default (dummy, 1, ' ')
          endif
       endif

    endif

    do m = 1,nwetdep

       call addfld (trim(wetdep_list(m))//'SFWET', horiz_only,  'A','kg/m2/s', &
            'Wet deposition flux at surface')
       call addfld (trim(wetdep_list(m))//'SFSIC', horiz_only,  'A','kg/m2/s', &
            'Wet deposition flux (incloud, convective) at surface')
       call addfld (trim(wetdep_list(m))//'SFSIS', horiz_only,  'A','kg/m2/s', &
            'Wet deposition flux (incloud, stratiform) at surface')
       call addfld (trim(wetdep_list(m))//'SFSBC', horiz_only,  'A','kg/m2/s', &
            'Wet deposition flux (belowcloud, convective) at surface')
       call addfld (trim(wetdep_list(m))//'SFSBS', horiz_only,  'A','kg/m2/s', &
            'Wet deposition flux (belowcloud, stratiform) at surface')
       call addfld (trim(wetdep_list(m))//'WET',   (/ 'lev' /), 'A','kg/kg/s', &
            'wet deposition tendency')
       call addfld (trim(wetdep_list(m))//'SIC',   (/ 'lev' /), 'A','kg/kg/s', &
            trim(wetdep_list(m))//' ic wet deposition')
       call addfld (trim(wetdep_list(m))//'SIS',   (/ 'lev' /), 'A','kg/kg/s', &
            trim(wetdep_list(m))//' is wet deposition')
       call addfld (trim(wetdep_list(m))//'SBC',   (/ 'lev' /), 'A','kg/kg/s', &
            trim(wetdep_list(m))//' bc wet deposition')
       call addfld (trim(wetdep_list(m))//'SBS',   (/ 'lev' /), 'A','kg/kg/s', &
            trim(wetdep_list(m))//' bs wet deposition')
    enddo

    if (nwetdep>0) then
       if (sslt_active) then
          dummy = 'SSTSFWET'
          call addfld (dummy,horiz_only, 'A','kg/m2/s','Sea salt wet deposition flux at surface')
          if ( history_aerosol ) then
             call add_default (dummy, 1, ' ')
          endif
       endif
       if (dust_active) then
          dummy = 'DSTSFWET'
          call addfld (dummy,horiz_only, 'A','kg/m2/s','Dust wet deposition flux at surface')
          if ( history_aerosol ) then
             call add_default (dummy, 1, ' ')
          endif
       endif
    endif

    if (dust_active) then
       ! emissions diagnostics ....

       do m = 1, dust_nbin
          dummy = trim(dust_names(m)) // 'SF'
          call addfld (dummy,horiz_only, 'A','kg/m2/s',trim(dust_names(m))//' dust surface emission')
          if (history_aerosol) then
             call add_default (dummy, 1, ' ')
          endif
       enddo

       dummy = 'DSTSFMBL'
       call addfld (dummy,horiz_only, 'A','kg/m2/s','Mobilization flux at surface')
       if (history_aerosol .or. history_dust) then
          call add_default (dummy, 1, ' ')
       endif

       dummy = 'LND_MBL'
       call addfld (dummy,horiz_only, 'A','frac','Soil erodibility factor')
       if (history_aerosol) then
          call add_default (dummy, 1, ' ')
       endif

    endif

    if (sslt_active) then

       dummy = 'SSTSFMBL'
       call addfld (dummy,horiz_only, 'A','kg/m2/s','Mobilization flux at surface')
       if (history_aerosol) then
          call add_default (dummy, 1, ' ')
       endif

       do m = 1, seasalt_nbin
          dummy = trim(seasalt_names(m)) // 'SF'
          call addfld (dummy,horiz_only, 'A','kg/m2/s',trim(seasalt_names(m))//' seasalt surface emission')
          if (history_aerosol) then
             call add_default (dummy, 1, ' ')
          endif
       enddo

    endif

    if( has_sox ) then
       call addfld( 'XPH_LWC',(/ 'lev' /), 'A','kg/kg', 'pH value multiplied by lwc')

       if ( history_aerosol ) then
          call add_default ('XPH_LWC', 1, ' ')
       endif
    endif

    so4_ndx    = get_spc_ndx( 'SO4' )
    soa_ndx    = get_spc_ndx( 'SOA' )
    soai_ndx   = get_spc_ndx( 'SOAI' )
    soam_ndx   = get_spc_ndx( 'SOAM' )
    soab_ndx   = get_spc_ndx( 'SOAB' )
    soat_ndx   = get_spc_ndx( 'SOAT' )
    soax_ndx   = get_spc_ndx( 'SOAX' )
    cb2_ndx    = get_spc_ndx( 'CB2' )
    oc2_ndx    = get_spc_ndx( 'OC2' )
    nit_ndx    = get_spc_ndx( 'NH4NO3' )

  end subroutine aero_model_init

  !=============================================================================
  !=============================================================================
  subroutine aero_model_drydep  ( state, pbuf, obklen, ustar, cam_in, dt, cam_out, ptend )

    use dust_sediment_mod, only: dust_sediment_tend
    use drydep_mod,        only: d3ddflux, calcram
    use dust_model,        only: dust_depvel, dust_nbin, dust_names
    use seasalt_model,     only: sslt_depvel=>seasalt_depvel, sslt_nbin=>seasalt_nbin, sslt_names=>seasalt_names

    ! args
    type(physics_state),    intent(in)    :: state     ! Physics state variables
    real(r8),               intent(in)    :: obklen(:)
    real(r8),               intent(in)    :: ustar(:)  ! sfc fric vel
    type(cam_in_t), target, intent(in)    :: cam_in    ! import state
    real(r8),               intent(in)    :: dt        ! time step
    type(cam_out_t),        intent(inout) :: cam_out   ! export state
    type(physics_ptend),    intent(out)   :: ptend     ! indivdual parameterization tendencies
    type(physics_buffer_desc),    pointer :: pbuf(:)

  ! local vars
    real(r8), pointer :: landfrac(:) ! land fraction
    real(r8), pointer :: icefrac(:)  ! ice fraction
    real(r8), pointer :: ocnfrac(:)  ! ocean fraction
    real(r8), pointer :: fvin(:)     !
    real(r8), pointer :: ram1in(:)   ! for dry dep velocities from land model for progseasalts

    real(r8) :: fv(pcols)            ! for dry dep velocities, from land modified over ocean & ice
    real(r8) :: ram1(pcols)          ! for dry dep velocities, from land modified over ocean & ice

     ! local decarations

    integer, parameter :: naero = sslt_nbin+dust_nbin
    integer, parameter :: begslt = 1
    integer, parameter :: endslt = sslt_nbin
    integer, parameter :: begdst = sslt_nbin+1
    integer, parameter :: enddst = sslt_nbin+dust_nbin

    integer :: ncol, lchnk

    character(len=6) :: aeronames(naero) ! = (/ sslt_names, dust_names /)

    real(r8) :: vlc_trb(pcols,naero)    !Turbulent deposn velocity (m/s)
    real(r8) :: vlc_grv(pcols,pver,naero)  !grav deposn velocity (m/s)
    real(r8) :: vlc_dry(pcols,pver,naero)  !dry deposn velocity (m/s)

    real(r8) :: dep_trb(pcols)       !kg/m2/s
    real(r8) :: dep_grv(pcols)       !kg/m2/s (total of grav and trb)

    real(r8) :: tsflx_dst(pcols)
    real(r8) :: tsflx_slt(pcols)
    real(r8) :: pvaeros(pcols,pverp)    ! sedimentation velocity in Pa
    real(r8) :: sflx(pcols)

    real(r8) :: tvs(pcols,pver)
    real(r8) :: rho(pcols,pver)      ! air density in kg/m3

    integer :: m,mm, i, im

    if (ndrydep<1) return

    landfrac => cam_in%landfrac(:)
    icefrac  => cam_in%icefrac(:)
    ocnfrac  => cam_in%ocnfrac(:)
    fvin     => cam_in%fv(:)
    ram1in   => cam_in%ram1(:)

    lchnk = state%lchnk
    ncol  = state%ncol

    ! calc ram and fv over ocean and sea ice ...
    call calcram( ncol,landfrac,icefrac,ocnfrac,obklen,&
                  ustar,ram1in,ram1,state%t(:,pver),state%pmid(:,pver),&
                  state%pdel(:,pver),fvin,fv)

    call outfld( 'airFV', fv(:), pcols, lchnk )
    call outfld( 'RAM1', ram1(:), pcols, lchnk )

    ! note that tendencies are not only in sfc layer (because of sedimentation)
    ! and that ptend is updated within each subroutine for different species

    call physics_ptend_init(ptend, state%psetcols, 'aero_model_drydep', lq=drydep_lq)

    aeronames(:sslt_nbin)   = sslt_names(:)
    aeronames(sslt_nbin+1:) = dust_names(:)

    lchnk = state%lchnk
    ncol  = state%ncol

    tvs(:ncol,:) = state%t(:ncol,:)
    rho(:ncol,:) = state%pmid(:ncol,:)/(rair*state%t(:ncol,:))

    ! compute dep velocities for sea salt and dust...
    if (sslt_active) then
       call sslt_depvel( state%t(:,:), state%pmid(:,:), state%q(:,:,1), ram1, fv, ncol, lchnk, &
                         vlc_dry(:,:,begslt:endslt), vlc_trb(:,begslt:endslt), vlc_grv(:,:,begslt:endslt))
    endif
    if (dust_active) then
       call dust_depvel( state%t(:,:), state%pmid(:,:),                 ram1, fv, ncol, &
                         vlc_dry(:,:,begdst:enddst), vlc_trb(:,begdst:enddst), vlc_grv(:,:,begdst:enddst) )
    endif

    tsflx_dst(:)=0._r8
    tsflx_slt(:)=0._r8

    ! do drydep for each of the bins of dust and seasalt
    do m=1,ndrydep

       mm = drydep_indices(m)
       findindex: do im = 1,naero
         if (trim(cnst_name(mm))==trim(aeronames(im))) exit findindex
       enddo findindex

       pvaeros(:ncol,1)=0._r8
       pvaeros(:ncol,2:pverp) = vlc_dry(:ncol,:,im)

       call outfld( trim(cnst_name(mm))//'DV', pvaeros(:,2:pverp), pcols, lchnk )

       if(.true.) then ! use phil's method
          !      convert from meters/sec to pascals/sec
          !      pvaeros(:,1) is assumed zero, use density from layer above in conversion
          pvaeros(:ncol,2:pverp) = pvaeros(:ncol,2:pverp) * rho(:ncol,:)*gravit

          !      calculate the tendencies and sfc fluxes from the above velocities
          call dust_sediment_tend( &
               ncol,             dt,       state%pint(:,:), state%pmid, state%pdel, state%t , &
               state%q(:,:,mm) , pvaeros  , ptend%q(:,:,mm), sflx  )
       else   !use charlie's method
          call d3ddflux(ncol, vlc_dry(:,:,im), state%q(:,:,mm),state%pmid,state%pdel, tvs,sflx,ptend%q(:,:,mm),dt)
       endif
       ! apportion dry deposition into turb and gravitational settling for tapes
       do i=1,ncol
          dep_trb(i)=sflx(i)*vlc_trb(i,im)/vlc_dry(i,pver,im)
          dep_grv(i)=sflx(i)*vlc_grv(i,pver,im)/vlc_dry(i,pver,im)
       enddo

       if ( any( sslt_names(:)==trim(cnst_name(mm)) ) ) &
            tsflx_slt(:ncol)=tsflx_slt(:ncol)+sflx(:ncol)
       if ( any( dust_names(:)==trim(cnst_name(mm)) ) ) &
            tsflx_dst(:ncol)=tsflx_dst(:ncol)+sflx(:ncol)

       ! if the user has specified prescribed aerosol dep fluxes then
       ! do not set cam_out dep fluxes according to the prognostic aerosols
       if (.not. aerodep_flx_prescribed()) then
          ! set deposition in export state
          if (im==begdst) then
             cam_out%dstdry1(:ncol) = max(sflx(:ncol), 0._r8)
          elseif(im==begdst+1) then
             cam_out%dstdry2(:ncol) = max(sflx(:ncol), 0._r8)
          elseif(im==begdst+2) then
             cam_out%dstdry3(:ncol) = max(sflx(:ncol), 0._r8)
          elseif(im==begdst+3) then
             cam_out%dstdry4(:ncol) = max(sflx(:ncol), 0._r8)
          endif
       endif

       call outfld( trim(cnst_name(mm))//'DD', sflx, pcols, lchnk)
       call outfld( trim(cnst_name(mm))//'TB', dep_trb, pcols, lchnk )
       call outfld( trim(cnst_name(mm))//'GV', dep_grv, pcols, lchnk )
       call outfld( trim(cnst_name(mm))//'DT', ptend%q(:,:,mm), pcols, lchnk)

    end do

    ! output the total dry deposition
    if (sslt_active) then
       call outfld( 'SSTSFDRY', tsflx_slt, pcols, lchnk)
    endif
    if (dust_active) then
       call outfld( 'DSTSFDRY', tsflx_dst, pcols, lchnk)
    endif

  endsubroutine aero_model_drydep

  !=============================================================================
  !=============================================================================
  subroutine aero_model_wetdep( state, dt, dlf, cam_out, ptend, pbuf)

    use wetdep,        only : wetdepa_v1, wetdep_inputs_t
    use dust_model,    only : dust_names
    use seasalt_model, only : sslt_names=>seasalt_names

    ! args

    type(physics_state), intent(in)    :: state       ! Physics state variables
    real(r8),            intent(in)    :: dt          ! time step
    real(r8),            intent(in)    :: dlf(:,:)    ! shallow+deep convective detrainment [kg/kg/s]
    type(cam_out_t),     intent(inout) :: cam_out     ! export state
    type(physics_ptend), intent(out)   :: ptend       ! indivdual parameterization tendencies
    type(physics_buffer_desc), pointer :: pbuf(:)

    ! local vars

    integer  :: ncol                     ! number of atmospheric columns
    integer  :: lchnk                    ! chunk identifier
    integer  :: m,mm, i,k

    real(r8) :: sflx_tot_dst(pcols)
    real(r8) :: sflx_tot_slt(pcols)

    real(r8) :: iscavt(pcols, pver)
    real(r8) :: scavt(pcols, pver)
    real(r8) :: scavcoef(pcols,pver)     ! Dana and Hales coefficient (/mm) (0.1)
    real(r8) :: sflx(pcols)              ! deposition flux

    real(r8) :: icscavt(pcols, pver)
    real(r8) :: isscavt(pcols, pver)
    real(r8) :: bcscavt(pcols, pver)
    real(r8) :: bsscavt(pcols, pver)

    real(r8) :: sol_factb, sol_facti

    real(r8) :: rainmr(pcols,pver)       ! mixing ratio of rain within cloud volume
    real(r8) :: cldv(pcols,pver)         ! cloudy volume undergoing scavenging
    real(r8) :: cldvcu(pcols,pver)       ! Convective precipitation area at the top interface of current layer
    real(r8) :: cldvst(pcols,pver)       ! Stratiform precipitation area at the top interface of current layer

    real(r8), pointer :: fracis(:,:,:)   ! fraction of transported species that are insoluble

    type(wetdep_inputs_t) dep_inputs  ! obj that contains inputs to wetdepa routine

    if (nwetdep<1) return

    call pbuf_get_field(pbuf, fracis_idx, fracis, start=(/1,1,1/), kount=(/pcols, pver, pcnst/) )

    call physics_ptend_init(ptend, state%psetcols, 'aero_model_wetdep', lq=wetdep_lq)

    call dep_inputs%init()
    call dep_inputs%set(state, pbuf)

    lchnk = state%lchnk
    ncol  = state%ncol

    sflx_tot_dst(:) = 0._r8
    sflx_tot_slt(:) = 0._r8

    do m = 1, nwetdep

       mm = wetdep_indices(m)

       sol_factb = aer_sol_factb(m)
       sol_facti = aer_sol_facti(m)

       scavcoef(:ncol,:) = aer_scav_coef(m)

       call wetdepa_v1( state%t, state%pmid, state%q(:,:,1), state%pdel, &
            dep_inputs%cldt, dep_inputs%cldcu, dep_inputs%cmfdqr, &
            dep_inputs%conicw, dep_inputs%prain, dep_inputs%qme, &
            dep_inputs%evapr, dep_inputs%totcond, state%q(:,:,mm), dt, &
            scavt, iscavt, dep_inputs%cldv, &
            fracis(:,:,mm), sol_factb, ncol, &
            scavcoef, &
            sol_facti_in=sol_facti, &
            icscavt=icscavt, isscavt=isscavt, bcscavt=bcscavt, bsscavt=bsscavt )

       ptend%q(:ncol,:,mm)=scavt(:ncol,:)

       call outfld( trim(cnst_name(mm))//'WET', ptend%q(:,:,mm), pcols, lchnk)
       call outfld( trim(cnst_name(mm))//'SIC', icscavt , pcols, lchnk)
       call outfld( trim(cnst_name(mm))//'SIS', isscavt, pcols, lchnk)
       call outfld( trim(cnst_name(mm))//'SBC', bcscavt, pcols, lchnk)
       call outfld( trim(cnst_name(mm))//'SBS', bsscavt, pcols, lchnk)

       sflx(:)=0._r8

       do k=1,pver
          do i=1,ncol
             sflx(i)=sflx(i)+ptend%q(i,k,mm)*state%pdel(i,k)/gravit
          enddo
       enddo
       call outfld( trim(cnst_name(mm))//'SFWET', sflx, pcols, lchnk)

       if ( any( sslt_names(:)==trim(cnst_name(mm)) ) ) &
            sflx_tot_slt(:ncol) = sflx_tot_slt(:ncol) + sflx(:ncol)
       if ( any( dust_names(:)==trim(cnst_name(mm)) ) ) &
            sflx_tot_dst(:ncol) = sflx_tot_dst(:ncol) + sflx(:ncol)

       ! if the user has specified prescribed aerosol dep fluxes then
       ! do not set cam_out dep fluxes according to the prognostic aerosols
       if (.not.aerodep_flx_prescribed()) then
          ! export deposition fluxes to coupler ??? why "-" sign ???
          if (trim(cnst_name(mm))=='CB2') then
             cam_out%bcphiwet(:) = max(-sflx(:), 0._r8)
          elseif (trim(cnst_name(mm))=='OC2') then
             cam_out%ocphiwet(:) = max(-sflx(:), 0._r8)
          elseif (trim(cnst_name(mm))==trim(dust_names(1))) then
             cam_out%dstwet1(:) = max(-sflx(:), 0._r8)
          elseif (trim(cnst_name(mm))==trim(dust_names(2))) then
             cam_out%dstwet2(:) = max(-sflx(:), 0._r8)
          elseif (trim(cnst_name(mm))==trim(dust_names(3))) then
             cam_out%dstwet3(:) = max(-sflx(:), 0._r8)
          elseif (trim(cnst_name(mm))==trim(dust_names(4))) then
             cam_out%dstwet4(:) = max(-sflx(:), 0._r8)
          endif
       endif

    enddo

    if (sslt_active) then
       call outfld( 'SSTSFWET', sflx_tot_slt, pcols, lchnk)
    endif
    if (dust_active) then
       call outfld( 'DSTSFWET', sflx_tot_dst, pcols, lchnk)
    endif

  endsubroutine aero_model_wetdep

  !-------------------------------------------------------------------------
  ! provides aerosol surface area info for sectional aerosols
  ! called from mo_usrrxt
  !-------------------------------------------------------------------------
  subroutine aero_model_surfarea( &
                  mmr, radmean, relhum, pmid, temp, strato_sad, sulfate,  m, ltrop, &
                  dlat, het1_ndx, pbuf, ncol, sfc, dm_aer, sad_total, reff_trop )

    use mo_constants, only : pi, avo => avogadro

    ! dummy args
    real(r8), intent(in)    :: pmid(:,:)
    real(r8), intent(in)    :: temp(:,:)
    real(r8), intent(in)    :: mmr(:,:,:)
    real(r8), intent(in)    :: radmean      ! mean radii in cm
    real(r8), intent(in)    :: strato_sad(:,:)
    integer,  intent(in)    :: ncol
    integer,  intent(in)    :: ltrop(:)
    real(r8), intent(in)    :: dlat(:)                    ! degrees latitude
    integer,  intent(in)    :: het1_ndx
    real(r8), intent(in)    :: relhum(:,:)
    real(r8), intent(in)    :: m(:,:) ! total atm density (/cm^3)
    real(r8), intent(in)    :: sulfate(:,:)
    type(physics_buffer_desc), pointer :: pbuf(:)

    real(r8), intent(inout) :: sfc(:,:,:)
    real(r8), intent(inout) :: dm_aer(:,:,:)
    real(r8), intent(inout) :: sad_total(:,:)
    real(r8), intent(out)   :: reff_trop(:,:)

    ! local vars

    integer  :: i,k
    real(r8) :: rho_air
    real(r8) :: v, n, n_exp, r_rd, r_sd
    real(r8) :: dm_sulf, dm_sulf_wet, log_sd_sulf, sfc_sulf, sfc_nit
    real(r8) :: dm_orgc, dm_orgc_wet, log_sd_orgc, sfc_oc, sfc_soa
    real(r8) :: sfc_soai, sfc_soam, sfc_soab, sfc_soat, sfc_soax
    real(r8) :: dm_bc, dm_bc_wet, log_sd_bc, sfc_bc
    real(r8) :: rxt_sulf, rxt_nit, rxt_oc, rxt_soa
    real(r8) :: c_n2o5, c_ho2, c_no2, c_no3
    real(r8) :: s_exp

    !-----------------------------------------------------------------
    ! 	... parameters for log-normal distribution by number
    ! references:
    !   Chin et al., JAS, 59, 461, 2003
    !   Liao et al., JGR, 108(D1), 4001, 2003
    !   Martin et al., JGR, 108(D3), 4097, 2003
    !-----------------------------------------------------------------
    real(r8), parameter :: rm_sulf  = 6.95e-6_r8        ! mean radius of sulfate particles (cm) (Chin)
    real(r8), parameter :: sd_sulf  = 2.03_r8           ! standard deviation of radius for sulfate (Chin)
    real(r8), parameter :: rho_sulf = 1.7e3_r8          ! density of sulfate aerosols (kg/m3) (Chin)

    real(r8), parameter :: rm_orgc  = 2.12e-6_r8        ! mean radius of organic carbon particles (cm) (Chin)
    real(r8), parameter :: sd_orgc  = 2.20_r8           ! standard deviation of radius for OC (Chin)
    real(r8), parameter :: rho_orgc = 1.8e3_r8          ! density of OC aerosols (kg/m3) (Chin)

    real(r8), parameter :: rm_bc    = 1.18e-6_r8        ! mean radius of soot/BC particles (cm) (Chin)
    real(r8), parameter :: sd_bc    = 2.00_r8           ! standard deviation of radius for BC (Chin)
    real(r8), parameter :: rho_bc   = 1.0e3_r8          ! density of BC aerosols (kg/m3) (Chin)

    real(r8), parameter :: mw_so4 = 98.e-3_r8     ! so4 molecular wt (kg/mole)

    integer  ::  irh, rh_l, rh_u
    real(r8) ::  factor, rfac_sulf, rfac_oc, rfac_bc, rfac_ss
    logical :: zero_aerosols

    !-----------------------------------------------------------------
    ! 	... table for hygroscopic growth effect on radius (Chin et al)
    !           (no growth effect for mineral dust)
    !-----------------------------------------------------------------
    real(r8), dimension(7) :: table_rh, table_rfac_sulf, table_rfac_bc, table_rfac_oc, table_rfac_ss

    data table_rh(1:7)        / 0.0_r8, 0.5_r8, 0.7_r8, 0.8_r8, 0.9_r8, 0.95_r8, 0.99_r8/
    data table_rfac_sulf(1:7) / 1.0_r8, 1.4_r8, 1.5_r8, 1.6_r8, 1.8_r8, 1.9_r8,  2.2_r8/
    data table_rfac_oc(1:7)   / 1.0_r8, 1.2_r8, 1.4_r8, 1.5_r8, 1.6_r8, 1.8_r8,  2.2_r8/
    data table_rfac_bc(1:7)   / 1.0_r8, 1.0_r8, 1.0_r8, 1.2_r8, 1.4_r8, 1.5_r8,  1.9_r8/
    data table_rfac_ss(1:7)   / 1.0_r8, 1.6_r8, 1.8_r8, 2.0_r8, 2.4_r8, 2.9_r8,  4.8_r8/

    !-----------------------------------------------------------------
    ! 	... exponent for calculating number density
    !-----------------------------------------------------------------
    n_exp = exp( -4.5_r8*log(sd_sulf)*log(sd_sulf) )

    dm_sulf = 2._r8 * rm_sulf
    dm_orgc = 2._r8 * rm_orgc
    dm_bc   = 2._r8 * rm_bc

    log_sd_sulf = log(sd_sulf)
    log_sd_orgc = log(sd_orgc)
    log_sd_bc   = log(sd_bc)

    reff_trop(:,:) = 0._r8

    ver_loop: do k = 1,pver
       col_loop: do i = 1,ncol
          !-------------------------------------------------------------------------
          ! 	... air density (kg/m3)
          !-------------------------------------------------------------------------
          rho_air = pmid(i,k)/(temp(i,k)*287.04_r8)
          !-------------------------------------------------------------------------
          !       ... aerosol growth interpolated from M.Chin's table
          !-------------------------------------------------------------------------
          if (relhum(i,k) >= table_rh(7)) then
             rfac_sulf = table_rfac_sulf(7)
             rfac_oc = table_rfac_oc(7)
             rfac_bc = table_rfac_bc(7)
          else
             do irh = 2,7
                if (relhum(i,k) <= table_rh(irh)) then
                   exit
                end if
             end do
             rh_l = irh-1
             rh_u = irh

             factor = (relhum(i,k) - table_rh(rh_l))/(table_rh(rh_u) - table_rh(rh_l))

             rfac_sulf = table_rfac_sulf(rh_l) + factor*(table_rfac_sulf(rh_u) - table_rfac_sulf(rh_l))
             rfac_oc = table_rfac_oc(rh_u) + factor*(table_rfac_oc(rh_u) - table_rfac_oc(rh_l))
             rfac_bc = table_rfac_bc(rh_u) + factor*(table_rfac_bc(rh_u) - table_rfac_bc(rh_l))
          end if

          dm_sulf_wet = dm_sulf * rfac_sulf
          dm_orgc_wet = dm_orgc * rfac_oc
          dm_bc_wet = dm_bc * rfac_bc

          dm_bc_wet   = min(dm_bc_wet  ,50.e-6_r8) ! maximum size is 0.5 micron (Chin)
          dm_orgc_wet = min(dm_orgc_wet,50.e-6_r8) ! maximum size is 0.5 micron (Chin)


          !-------------------------------------------------------------------------
          ! 	... sulfate aerosols
          !-------------------------------------------------------------------------
          zero_aerosols = k < ltrop(i)
          if ( abs( dlat(i) ) > 50._r8 ) then
             zero_aerosols = pmid(i,k) < 30000._r8
          endif
          !-------------------------------------------------------------------------
          !       ... use ubvals climatology for stratospheric sulfate surface area density
          !-------------------------------------------------------------------------
          if( zero_aerosols ) then
             sfc_sulf = strato_sad(i,k)
             if ( het1_ndx > 0 ) then
                sfc_sulf = 0._r8        ! reaction already taken into account in mo_strato_rates.F90
             end if
             sfc_nit = 0._r8
             sfc_soa = 0._r8
             sfc_oc  = 0._r8
             sfc_bc  = 0._r8
          else

             if( so4_ndx > 0 ) then
                !-------------------------------------------------------------------------
                ! convert mass mixing ratio of aerosol to cm3/cm3 (cm^3_aerosol/cm^3_air)
                ! v=volume density (m^3/m^3)
                ! rho_aer=density of aerosol (kg/m^3)
                ! v=m*rho_air/rho_aer   [kg/kg * (kg/m3)_air/(kg/m3)_aer]
                !-------------------------------------------------------------------------
                v = mmr(i,k,so4_ndx) * rho_air/rho_sulf
                !-------------------------------------------------------------------------
                ! calculate the number density of aerosol (aerosols/cm3)
                ! assuming a lognormal distribution
                ! n  = (aerosols/cm3)
                ! dm = geometric mean diameter
                !
                ! because only the dry mass of the aerosols is known, we
                ! use the mean dry radius
                !-------------------------------------------------------------------------
                n  = v * (6._r8/pi)*(1._r8/(dm_sulf**3._r8))*n_exp
                !-------------------------------------------------------------------------
                ! find surface area of aerosols using dm_wet, log_sd
                !  (increase of sd due to RH is negligible)
                ! and number density calculated above as distribution
                ! parameters
                ! sfc = surface area of wet aerosols (cm^2/cm^3)
                !-------------------------------------------------------------------------
                s_exp    = exp(2._r8*log_sd_sulf*log_sd_sulf)
                sfc_sulf = n * pi * (dm_sulf_wet**2._r8) * s_exp

             else
                !-------------------------------------------------------------------------
                !  if so4 not simulated, use off-line sulfate and calculate as above
                !  convert sulfate vmr to volume density of aerosol (cm^3_aerosol/cm^3_air)
                !-------------------------------------------------------------------------
                v = sulfate(i,k) * m(i,k) * mw_so4 / (avo * rho_sulf) *1.e6_r8
                n  = v * (6._r8/pi)*(1._r8/(dm_sulf**3._r8))*n_exp
                s_exp    = exp(2._r8*log_sd_sulf*log_sd_sulf)
                sfc_sulf = n * pi * (dm_sulf_wet**2._r8) * s_exp

             end if

             !-------------------------------------------------------------------------
             ! ammonium nitrate (follow same procedure as sulfate, using size and density of sulfate)
             !-------------------------------------------------------------------------
             if( nit_ndx > 0 ) then
                v = mmr(i,k,nit_ndx) * rho_air/rho_sulf
                n  = v * (6._r8/pi)*(1._r8/(dm_sulf**3._r8))*n_exp
                s_exp   = exp(2._r8*log_sd_sulf*log_sd_sulf)
                sfc_nit = n * pi * (dm_sulf_wet**2._r8) * s_exp
             else
                sfc_nit = 0._r8
             end if

             !-------------------------------------------------------------------------
             ! hydrophylic organic carbon (follow same procedure as sulfate)
             !-------------------------------------------------------------------------
             if( oc2_ndx > 0 ) then
                v = mmr(i,k,oc2_ndx) * rho_air/rho_orgc
                n  = v * (6._r8/pi)*(1._r8/(dm_orgc**3))*n_exp
                s_exp    = exp(2._r8*log_sd_orgc*log_sd_orgc)
                sfc_oc   = n * pi * (dm_orgc_wet**2._r8) * s_exp
             else
                sfc_oc = 0._r8
             end if

             !-------------------------------------------------------------------------
             ! secondary organic carbon (follow same procedure as sulfate)
             !-------------------------------------------------------------------------
             if( soa_ndx > 0 ) then
                v = mmr(i,k,soa_ndx) * rho_air/rho_orgc
                n  = v * (6._r8/pi)*(1._r8/(dm_orgc**3._r8))*n_exp
                s_exp     = exp(2._r8*log_sd_orgc*log_sd_orgc)
                sfc_soa   = n * pi * (dm_orgc_wet**2._r8) * s_exp
             else
                sfc_soa = 0._r8
             end if

             !-------------------------------------------------------------------------
             ! black carbon (follow same procedure as sulfate)
             !-------------------------------------------------------------------------
             if( cb2_ndx > 0 ) then
                v = mmr(i,k,cb2_ndx) * rho_air/rho_bc
                n  = v * (6._r8/pi)*(1._r8/(dm_bc**3._r8))*n_exp
                s_exp     = exp(2._r8*log_sd_bc*log_sd_bc)
                sfc_bc   = n * pi * (dm_bc_wet**2._r8) * s_exp
             else
                sfc_bc = 0._r8
             end if
             if( soai_ndx > 0 ) then
                v = mmr(i,k,soai_ndx) * rho_air/rho_orgc
                n  = v * (6._r8/pi)*(1._r8/(dm_orgc**3._r8))*n_exp
                s_exp     = exp(2._r8*log_sd_orgc*log_sd_orgc)
                sfc_soai   = n * pi * (dm_orgc_wet**2._r8) * s_exp
             else
                sfc_soai = 0._r8
             end if
             if( soam_ndx > 0 ) then
                v = mmr(i,k,soam_ndx) * rho_air/rho_orgc
                n  = v * (6._r8/pi)*(1._r8/(dm_orgc**3._r8))*n_exp
                s_exp     = exp(2._r8*log_sd_orgc*log_sd_orgc)
                sfc_soam   = n * pi * (dm_orgc_wet**2._r8) * s_exp
             else
                sfc_soam = 0._r8
             end if
             if( soab_ndx > 0 ) then
                v = mmr(i,k,soab_ndx) * rho_air/rho_orgc
                n  = v * (6._r8/pi)*(1._r8/(dm_orgc**3._r8))*n_exp
                s_exp     = exp(2._r8*log_sd_orgc*log_sd_orgc)
                sfc_soab   = n * pi * (dm_orgc_wet**2._r8) * s_exp
             else
                sfc_soab = 0._r8
             end if
             if( soat_ndx > 0 ) then
                v = mmr(i,k,soat_ndx) * rho_air/rho_orgc
                n  = v * (6._r8/pi)*(1._r8/(dm_orgc**3._r8))*n_exp
                s_exp     = exp(2._r8*log_sd_orgc*log_sd_orgc)
                sfc_soat   = n * pi * (dm_orgc_wet**2._r8) * s_exp
             else
                sfc_soat = 0._r8
             end if
             if( soax_ndx > 0 ) then
                v = mmr(i,k,soax_ndx) * rho_air/rho_orgc
                n  = v * (6._r8/pi)*(1._r8/(dm_orgc**3._r8))*n_exp
                s_exp     = exp(2._r8*log_sd_orgc*log_sd_orgc)
                sfc_soax   = n * pi * (dm_orgc_wet**2._r8) * s_exp
             else
                sfc_soax = 0._r8
             end if
             sfc_soa = sfc_soa + sfc_soai + sfc_soam + sfc_soab + sfc_soat + sfc_soax

          end if

          sfc(i,k,:) = (/ sfc_sulf, sfc_nit, sfc_oc, sfc_soa, sfc_bc /)
          dm_aer(i,k,:) = (/ dm_sulf_wet,dm_sulf_wet,dm_orgc_wet,dm_orgc_wet,dm_bc_wet /)

          !-------------------------------------------------------------------------
          !  	... add up total surface area density for output
          !-------------------------------------------------------------------------
          sad_total(i,k) = sfc_sulf + sfc_nit + sfc_oc + sfc_soa + sfc_bc

       enddo col_loop
    enddo ver_loop

  end subroutine aero_model_surfarea

  !-------------------------------------------------------------------------
  ! stub
  !-------------------------------------------------------------------------
  subroutine aero_model_strat_surfarea( ncol, mmr, pmid, temp, ltrop, pbuf, strato_sad, reff_strat )

    ! dummy args
    integer,  intent(in)    :: ncol
    real(r8), intent(in)    :: mmr(:,:,:)
    real(r8), intent(in)    :: pmid(:,:)
    real(r8), intent(in)    :: temp(:,:)
    integer,  intent(in)    :: ltrop(:) ! tropopause level indices
    type(physics_buffer_desc), pointer :: pbuf(:)
    real(r8), intent(out)   :: strato_sad(:,:)
    real(r8), intent(out)   :: reff_strat(:,:)

    strato_sad(:,:) = 0._r8
    reff_strat(:,:) = 0._r8

  end subroutine aero_model_strat_surfarea

  !=============================================================================
  !=============================================================================
  subroutine aero_model_gasaerexch( loffset, ncol, lchnk, troplev, delt, reaction_rates, &
                                    tfld, pmid, pdel, mbar, relhum, &
                                    zm,  qh2o, cwat, cldfr, cldnum, &
                                    airdens, invariants, del_h2so4_gasprod,  &
                                    vmr0, vmr, pbuf )

    use chem_mods,   only : gas_pcnst
    use mo_aerosols, only : aerosols_formation, has_aerosols
    use mo_setsox,   only : setsox, has_sox
    use mo_setsoa,   only : setsoa, has_soa

    !-----------------------------------------------------------------------
    !      ... dummy arguments
    !-----------------------------------------------------------------------
    integer,  intent(in) :: loffset                ! offset applied to modal aero "pointers"
    integer,  intent(in) :: ncol                   ! number columns in chunk
    integer,  intent(in) :: lchnk                  ! chunk index
    integer,  intent(in) :: troplev(:)
    real(r8), intent(in) :: delt                   ! time step size (sec)
    real(r8), intent(in) :: reaction_rates(:,:,:)  ! reaction rates
    real(r8), intent(in) :: tfld(:,:)              ! temperature (K)
    real(r8), intent(in) :: pmid(:,:)              ! pressure at model levels (Pa)
    real(r8), intent(in) :: pdel(:,:)              ! pressure thickness of levels (Pa)
    real(r8), intent(in) :: mbar(:,:)              ! mean wet atmospheric mass ( amu )
    real(r8), intent(in) :: relhum(:,:)            ! relative humidity
    real(r8), intent(in) :: airdens(:,:)           ! total atms density (molec/cm**3)
    real(r8), intent(in) :: invariants(:,:,:)
    real(r8), intent(in) :: del_h2so4_gasprod(:,:)
    real(r8), intent(in) :: zm(:,:)
    real(r8), intent(in) :: qh2o(:,:)
    real(r8), intent(in) :: cwat(:,:)          ! cloud liquid water content (kg/kg)
    real(r8), intent(in) :: cldfr(:,:)
    real(r8), intent(in) :: cldnum(:,:)       ! droplet number concentration (#/kg)
    real(r8), intent(in) :: vmr0(:,:,:)       ! initial mixing ratios (before gas-phase chem changes)
    real(r8), intent(inout) :: vmr(:,:,:)         ! mixing ratios ( vmr )

    type(physics_buffer_desc), pointer :: pbuf(:)

    ! local vars

    real(r8) :: vmrcw(ncol,pver,gas_pcnst)            ! cloud-borne aerosol (vmr)

    real(r8) ::  aqso4(ncol,1)               ! aqueous phase chemistry
    real(r8) ::  aqh2so4(ncol,1)             ! aqueous phase chemistry
    real(r8) ::  aqso4_h2o2(ncol)            ! SO4 aqueous phase chemistry due to H2O2
    real(r8) ::  aqso4_o3(ncol)              ! SO4 aqueous phase chemistry due to O3
    real(r8) ::  xphlwc(ncol,pver)           ! pH value multiplied by lwc


  ! aqueous chemistry ...

    if( has_sox ) then
       call setsox(   &
            ncol,     &
            lchnk,    &
            loffset,  &
            delt,     &
            pmid,     &
            pdel,     &
            tfld,     &
            mbar,     &
            cwat,     &
            cldfr,    &
            cldnum,   &
            airdens,  &
            invariants, &
            vmrcw,    &
            vmr,      &
            xphlwc,   &
            aqso4,    &
            aqh2so4,  &
            aqso4_h2o2,&
            aqso4_o3  &
            )
        call outfld( 'XPH_LWC',xphlwc(:ncol,:), ncol , lchnk )
    endif

    if( has_soa ) then
       call setsoa( ncol, lchnk, delt, reaction_rates, tfld, airdens, vmr, pbuf)
    endif

    if( has_aerosols ) then
       call aerosols_formation( ncol, lchnk, tfld, relhum, vmr )
    endif


  end subroutine aero_model_gasaerexch

  !=============================================================================
  !=============================================================================
  subroutine aero_model_emissions( state, cam_in )
    use seasalt_model, only: seasalt_emis, seasalt_indices
    use dust_model,    only: dust_emis, dust_indices
    use physics_types, only: physics_state

    ! Arguments:

    type(physics_state),    intent(in)    :: state   ! Physics state variables
    type(cam_in_t),         intent(inout) :: cam_in  ! import state

    ! local vars

    integer :: lchnk, ncol
    integer :: m, mm
    real(r8) :: soil_erod_tmp(pcols)
    real(r8) :: sflx(pcols)   ! accumulate over all bins for output
    real(r8) :: u10cubed(pcols)
    real (r8), parameter :: z0=0.0001_r8  ! m roughness length over oceans--from ocean model

    lchnk = state%lchnk
    ncol = state%ncol

    if (dust_active) then

       call dust_emis( ncol, lchnk, cam_in%dstflx, cam_in%cflx, soil_erod_tmp )

       ! some dust emis diagnostics ...
       sflx(:)=0._r8
       do m=1,dust_nbin
          mm = dust_indices(m)
          sflx(:ncol)=sflx(:ncol)+cam_in%cflx(:ncol,mm)
          call outfld(trim(dust_names(m))//'SF',cam_in%cflx(:,mm),pcols, lchnk)
       enddo
       call outfld('DSTSFMBL',sflx(:),pcols,lchnk)
       call outfld('LND_MBL',soil_erod_tmp(:),pcols, lchnk )
    endif

    if (sslt_active) then
       u10cubed(:ncol)=sqrt(state%u(:ncol,pver)**2+state%v(:ncol,pver)**2)
       ! move the winds to 10m high from the midpoint of the gridbox:
       ! follows Tie and Seinfeld and Pandis, p.859 with math.

       u10cubed(:ncol)=u10cubed(:ncol)*log(10._r8/z0)/log(state%zm(:ncol,pver)/z0)

       ! we need them to the 3.41 power, according to Gong et al., 1997:
       u10cubed(:ncol)=u10cubed(:ncol)**3.41_r8

       sflx(:)=0._r8

       call seasalt_emis( u10cubed, cam_in%sst, cam_in%ocnfrac, ncol, cam_in%cflx )

       do m=1,seasalt_nbin
          mm = seasalt_indices(m)
          sflx(:ncol)=sflx(:ncol)+cam_in%cflx(:ncol,mm)
          call outfld(trim(seasalt_names(m))//'SF',cam_in%cflx(:,mm),pcols,lchnk)
       enddo
       call outfld('SSTSFMBL',sflx(:),pcols,lchnk)
    endif

  end subroutine aero_model_emissions

end module aero_model
